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1  Introduction 

We  will  not  attempt  any  comprehensive  review  of  infinite  elements.  The  technology  has  long  been  in  place, 
has  enjoyed  a  tremendous  progress  and  a  large  number  of  research  papers  in  the  recent  years,  and  it  remains 
at  forefront  of  challenging  research,  see  e.g.  [2]. 

This  short  contribution  has  been  motivated  with  the  idea  of  the  second  author  who,  several  years  ago, 
suggested  to  employ  rational  functions  generated  from  the  (integrated)  Jacobi  polynomials  as  infinite  el¬ 
ement  shape  functions.  The  issue  resurfaced  recently  in  context  of  using  the  infinite  element  for  solving 
exterior  Maxwell  problems  and  determining  Radar  Cross  Sections  (RCS)  [13].  Apparently,  the  IE  elements 
have  suffered  from  poor  conditioning,  especially  for  non-spherical  truncating  surfaces.  The  issue  has  been 
reported  in  the  acoustics  community  as  well,  see  e.g.  [11]. 
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In  process  of  putting  together  a  short  chapter  on  /ip-adaptivity  for  exterior  boundary  value  problems  [4], 
the  first  author  has  learned  a  few  extra  new  (at  least  for  him)  facts  that  might  be  useful  to  the  IE  community. 
Those  include  the  following  issues. 

•  The  possibility  of  using  the  same  shape  trial  and  test  functions  for  the  conjugated  infinite  element. 
To  our  best  knowledge,  a  complete  convergence  proof  for  any  infinite  element  technology  remains  an 
open  issue  (cf.  [5]).  Symmetry  of  the  formulation  makes,  in  general,  the  stability  analysis  easier,  and 
can  be  exploited  in  designing  optimal  linear  solvers. 

•  The  possibility  of  a  uniform  formulation  for  both  2D  and  3D  problems. 

•  The  possibility  of  varying  locally  the  IE  spectral  order  of  approximation. 

•  The  possibility  of  a  direct  evaluation  of  RCS,  without  any  costly  post-processing. 

The  final  mofivafion  came  from  an  early  invesfigafion  of  fhe  effecl  of  resolving  singularities  in  fhe 
field  on  fhe  accuracy  of  RCS  calculations,  aimed  af  developing  a  general  fechnology  combining  fhe  goal- 
driven  automatic  /ip-adapfivify  [12,  10]  wifh  infinife  elemenfs.  The  aufomafic  /ip-adapfivify  allows  for  an 
unprecedenfed  accuracy  of  simulafions.  The  question  is  :  “Is  fhis  fechnology  really  necessary  for  accurafe 
RCS  predicfions?” 

2  Formulation  of  the  problem 

We  shall  consider  fhe  exferior  boundary  value  problem  for  fhe  Helmholtz  equation,  representing  a  class  of 
acoustical  rigid  scattering  problems. 

Eet  Uint  be  an  open  bounded  domain  inM^,  we  set  —  0,int-  Given  an  incident  pressure  wave 

ytnc  satisfies  fhe  Helmholfz  equafion  in  fhe  whole  space, 

=  0  mR^  , 

we  look  for  a  scattered  pressure  wave  u  fhaf  safisfies: 

•  fhe  same  Helmholtz  equation,  but  only  in  the  exterior  domain, 

—Art  —  k^u  =  0  in  =  .E”  —  , 

•  Neumann  boundary  condition  on  boundary  V  =  dkt  =  dklint, 

d{u  +  ^ 

•  Sommerfeld  radiation  condition  at  infinity, 

I 

—  +  iku  e  L^{n)  . 
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Here  A:  =  ^  >  0  is  the  wave  number  with  a;  representing  the  angular  frequeney,  and  c  the  speed  of 
sound  in  the  (homogeneous)  medium,  n  is  the  outward  unit  veetor,  and  r  stands  for  the  radial  eoordinate 
eorresponding  to  polar  eoordinates  with  the  origin  at  an  arbitrary  point,  typieally  ehosen  inside  of  the  interior 
domain  ^int  oeeupied  by  the  rigid  obstaele  (seatterer).  As  the  gradient  of  the  pressure  is  proportional  to  the 
veloeity  veetor,  the  Neumann  boundary  eondition  eorresponds  to  the  requirement  that  the  normal  eomponent 
of  the  total  veloeity  veetor  must  vanish  on  the  boundary  of  the  rigid  obstaele.  The  sum  of  the  ineident  and 
seattered  pressures  represents  the  total  pressure  wave.  In  praetiee,  the  ineident  pressure  wave  eorresponds 
to  a  souree  loeated  away  from  the  obstaele.  If  the  loeation  of  the  souree  is  far  enough,  the  ineident  pressure 
wave  may  be  assumed  in  the  form  of  a  plane  wave  (souree  at  infinity). 


where  uq  is  a  sealing  faetor,  and  unit  veetor  e  represents  the  direetion  of  the  ineoming  ineident  wave.  As 
the  problem  is  linear,  we  ean  assume  ttQ  =  1.  The  role  of  the  ineident  wave  is  therefore  to  provide  only  a 
driving  Neumann  load  data  on  boundary  T, 


du  _  _ 

dn  ^  dn 


(2.1) 


The  Sommerfeld  radiation  eondition  represents  a  requirement  that  the  seattered  wave  must  be  outgoing  to 
infinity,  and  it  ean  be  represented  in  many  different  but  equivalent  forms.  The  form  used  here  has  a  partieular 
advantage  of  being  independent  of  spaee  dimension. 

Notiee  that  the  sign  in  the  Sommerfeld  radiation  eondition  eorresponds  to  the  ansatz  in  time. 

Of  partieular  interest  to  us,  will  be  the  eomputation  of  the  monostatie  Radar  Cross  Section  (RCS)  defined 
as, 

lim  \u{rx)  |r“  .  (2.2) 

r^oo 

Here  ®  is  a  poinf  on  fhe  unif  sphere  (eirele)  speeifying  fhe  direefion  of  fhe  ineidenf  wave,  e  =  x,  and  u  is 
fhe  eorresponding  seaffered  wave.  Coeffieienf, 


a  = 


{ 


1 

2 

1 


refleefs  fhe  deeay  rafe  of  fhe  leading  term  (far  field  pattern)  in  fhe  solufion. 


3  Infinite  element  discretization 

3.1  Scattering  on  a  cylinder  (sphere) 

Infrodueing  fhe  variational  formulation  for  an  exterior  problem  is  differenf  from  fhe  sfandard  eonsfruefion 
for  bounded  domains.  We  will  begin  our  diseussion  wifh  a  speeial  ease  of  a  eylindrieal  (spherieal  in  3D) 
seatterer  of  radius  a. 
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We  surround  the  seatterer  with  a  larger  eirele  (sphere  in  3D)  of  radius  R,  and  follow  the  standard  deriva¬ 
tion  of  the  variational  formulation  for  the  truneated  exterior  domain, 


Or  =  {a:  :  a  <  |a:|  <  R}  . 


Assuming  u,v  E  we  obtain, 


VttVu 


k'^uv}  dx  -f  ik 


uv  ds 


L 


gv  ds  + 


wv  ds, 


Vu  . 


(3.3) 


Here  w  :=  ^  +  iku,  aeeording  to  the  Sommerfeld  eondition,  denotes  an  unknown  funetion,  and  T  is 
the  boundary  of  the  seatterer  i.e.,  in  this  ease,  T  =  Sa- 

We  want  to  pass  R  -foo  in  (3.3).  L^-integrability  of  funetion  w  implies  that  the  unknown  term  will 
vanish.  The  first  boundary  integral  on  the  right-hand  side  is  independent  of  R,  so  we  only  need  to  foeus  on 
the  two  integrals  on  the  left-hand  side  only. 


We  shall  use  the  following  eurvilinear  eoordinates,  eorresponding  to  the  finite  element  implementation. 


x{r,^)  =  rXaiO,  r>l,  |a:a|=a, 


where  Xa{^)  denotes  a  parametrization  of  the  eirele  (sphere  in  3D)  of  radius  a.  We  have  the  standard 
formulas, 

dx 

=  Xa  =  acr,  where  :=  a“^a:a, 
or 

dx  dXa  ,dXa,  ,  .dx.^dx 

=  r  e^,  where  :=  |  — | 


d^  d^  '  d^ 

_  1  du  1  du , dxn  ,  1 

Vu=  -—e, + 

a  dr  r  d^  d^ 


du 


dC  dC 


1. 


with  ^ Sa  denoting  the  gradient  on  the  eirele  (sphere)  Sa- 


Separation  of  variables 


We  shall  exploit  the  separability  of  the  domain  by  assuming  the  solution  u  and  the  test  funetion  v  in  the 
standard  form, 

u{r,Xa)  =  '^U,{r)'llj,{Xa),  v{r,Xa)  =  '^V^{r)(j)^{Xa)  . 

L  K 

In  the  following,  we  will  drop  indiees  i,  k  and  foeus  on  the  integrability  of  a  single  term.  On  the  diserete 
level,  eaeh  term  will  eorrespond  to  an  entry  in  the  global  stiffness  matrix. 

Separation  of  variables  leads  to  the  following  formula  for  the  domain  integral. 


a 


-1 


(kafUV)  r"dr  +  ika  p"  U(p)V(p)  }  /  dSa 


+a 


J  JSa 

r  \uVr^dr  [  Vsai^'^Sa^dSa, 
■h  ^  J  Sa 
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where  pa  =  R,  a  =  1  for  2D,  and  a  =  2  for  3D  problems.  We  shall  diseuss  now  suffieient  eonditions  for 
Lebesgue  integrability  of  terms, 

POO 


{u'V 


{kafUV)  r^dr, 


^UV  r^dr 


(3.4) 


and  the  existenee  of  the  limit. 


lim 

p—>oo 


U{p)V{p) 


(3.5) 


Step  1:  Taking  out  the  oscillating  factor.  The  formula  for  the  exaet  solution  (ef.  [9])  suggests  building 
into  the  solution  ansatz,  the  oseillating  faetor  “Overloading”  symbols  U  and  V,  we  introduee, 

U{r)  =  V{r)  = 

U'  =  {U'  -  V'  =  (V'  +  . 

This  leads  to  the  following  modiheations  of  the  integrands  (3.4)  and  the  boundary  term  (3.5), 

U'V'  -  {kafUV  =  U'V  +  ika{U'V  -  UV) 

Luv  =\uv 

/jp  ^ 

U{p)V{p)  =  U{p)V{p)  . 


step  2:  Taking  out  the  Jacobian.  In  an  effort  to  make  the  formulation  dimension  independent,  our  seeond 
modifieation  of  the  ansatz  for  the  solution  and  the  test  funetion,  involves  an  additional  faetor  eorresponding 
to  the  Jaeobian  r“,  where  a  =  1  for  2D,  and  a  =  2  for  3D  problems.  “Overloading”  symbols  U,  V  again, 
we  have, 

U  =  V  = 

U'  =  {U'  -  V  =  {V  -  — F)r-“/2  . 

This  leads  to  the  following  modifieation  of  the  terms  of  interest, 

2 

{u'V  +  ika{U'V  -  UV)}  r“  =  U'V  -  ^{UV)'  +  ^UV  +  ika{U'V  -  UV) 

p“U(p)V(p)  =U(p)V(p) 

Term  {UV)'  ean  be  eliminated  by  integrating  by  parts. 


a 

2  A  r 


—  rv  1  _  ty  _ 

{UVy  dr  =  --  ^UV  dr  +  -U{1)V{1)  . 
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The  final  form  of  the  two  terms  eontributing  to  the  domain  integral  beeomes, 


I"  U'V'  dr  +  1(1  -  l)  1  UV  dr  +  |t7(l)F(l) 

■  poc 

J  {U'V  -  UV')  dr  +  lim  U{p)V{p) 


+ika 


p—^oo 


(3.6) 


and, 


/ 


1 


1 


jUVdr, 


(3.7) 


where  symbols  U  and  V  eorrespond  to  the  eontributions  on  the  right  hand  side  of  the  formulas, 

U{r)  =  , 

V{r)  =  r-“/2e*^a(»-i)y(^)  . 


Step  3:  Selecting  radial  shape  functions.  The  formula  for  the  exaet  solution  suggests  seleeting, 

U  «  U]\f  =  tto  +  —ui  +  . . .  +  py  rtiv  ,  (3.8) 

where  rti, . . . ,  ttiv  are  unknown  eoeffieients. 

The  test  funetions  may  be  seleeted  in  different  ways.  We  shall  diseuss  two  most  natural  ehoiees. 


Same  trial  and  test  functions  (Bubnov-Galerkin) 

V  «  FiV  =  UO  +  -Vl  +  .  .  .  +  py  t’iV  5 

Notice  that, 

•  all  integrals  in  (3.6)  and  (3.7)  are  well  defined; 

•  \iinp^ooU{p)V{p)  =uqVq 

Different  trial  and  test  functions  (Petrov-Galerkin) 


In  this  case, 

•  all  integrals  in  (3.6)  and  (3.7)  are  well  defined; 

•  limp^oo(7(p)F(p)  =  0 

The  symmetric  formulation  is  possible  due  to  the  interpretation  of  the  domain  integral  in  the  Cauchy 
Principal  Value  (CPV)  sense.  We  first  integrate  over  the  truncated  domain,  then  allow  for  the  cancellation  of 
Lebesgue  non-summable  terms  (cf.  Step  1  of  the  derivation),  and  only  then  pass  to  the  limit  with  p  ^  oo.  For 
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the  second  choice  of  the  test  functions,  all  integrals  contributing  to  the  sesquilinear  form  exist  independently. 
In  the  three-  dimensional  case,  we  can  also  make  an  intermediate  choice, 

F  «  V/v  =  -  I  ^^0  +  +  •  •  •  + 

r  [  r 

In  this  case,  the  integrals  have  to  also  be  understood  in  the  CPV  sense. 

We  will  present  numerical  results  for  both  formulations. 


3.2  Conditioning  and  optimal  shape  functions. 


Selecting  functions  1  /r”  for  radial  shape  functions  results  in  poor  conditioning  as  the  functions  are  “too 
similar”  to  each  other.  In  order  to  make  a  better  choice,  we  first  make  the  change  of  variables, 

1  ,  1  ,  dC  ndU 

r  =  — ,  or  = - 2  dx,  =  —x  . 


X 


dr 


dx 


Then,  integrals  (3.6),  (3.7)  transform  into. 


L 


O', 


—ika 


f\u' 

Jo 


/  UV  dx+ -U{1)V{1) 
Iq  2 


V  -  UV)  dx  -  C(0)C(0) 


(3.9) 


and. 


L 


1 


UVdx  . 


(3.10) 


The  leading  term  (derivative  times  derivative)  is  most  responsible  for  the  conditioning  problems,  and  it 
motivates  our  choice  of  radial  shape  functions  to  be  integrals  of  Jacobi  polynomials  {^x  —  1),  which, 
thanks  to  the  orthogonality  of  the  Jacobi  polynomials,  are  mutually  orthogonal  with  respect  to  the  first 
integral  in  (3.9)  and  hence  will  lead  to  better-conditioned  linear  system.  We  note  that  the  approximation 
properties  of  the  associated  (through  the  transform  x  =  ^)  rational  functions  Rn{r)  :=  —  1)  and 

their  applications  for  solving  elliptic  equations  in  exterior  domains  were  studied  in  [8]. 


3.3  Scatterers  of  an  arbitrary  shape. 

For  a  scatterer  of  an  arbitrary  shape,  the  formulation  is  identical  to  the  special  scatter  considered  above. 
Meshing  the  exterior  domain  involves  introducing  a  truncating  sphere  (circle)  of  radius  a  -  Sa  with  a  usually 
selected  in  such  a  way  that  the  distance  of  Sa  from  the  scatterer  does  not  exceed  one  wave  length.  The  space 
between  the  scatterer  and  the  truncating  sphere  is  then  meshed  with  finite  element.  Each  finite  element  face 
(edge  in  2D)  on  the  truncating  sphere  gives  then  rise  to  a  separate  infinite  element.  The  technique  can  be 
extended  to  ellipsoidal  (elliptic)  truncating  surfaces,  see  the  illustration  of  the  concept  in  Fig.  1. 
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Figure  1 :  FE/IE  mesh  around  a  seatterer 


4  Implementation 

We  proeeed  now  with  a  short  deseription  of  the  implementation  of  the  deseribed  diseretization  within  the 
2Dhp  -  the  two-dimensional  hp  eode,  supporting  automatie  /ip-adaptivity  [7,  4]. 

4.1  Data  structure,  interface  with  frontal  solver 

Conceptually,  infinite  elements  may  be  treated  the  same  way  as  finite  elements.  Shape-wise,  in  2D,  they  fall 
into  the  category  of  quads.  The  corresponding  shape  functions  are  defined  as  fensor  producfs  of  ID  shape 
funcfions,  same  way  as  for  fhe  sfandard  quadrilaferals,  excepf  for  fhe  logic  on  nodes  and  fhe  corresponding 
nodal  connecfivifies.  In  2D,  each  infinite  elemenf  has  only  fwo  verfex  nodes,  one  mid-edge  node  on  fhe 
fruncafing  circle,  fhe  middle  node,  and  fwo  mid-edge  nodes  on  fhe  radial  edges  originafing  from  fhe  vertex 
nodes  on  fhe  fruncafing  circle. 

Thus,  one  way  fo  code  fhe  infinite  elemenfs,  would  be  fo  add  a  new  elemenf  fype  fo  fhe  code,  and 
develop  fhen  all  necessary  dafa  sfrucfure,  consfrained  approximafion  (in  3D  only),  mesh  refinemenfs,  and 
mesh  opfimizafion  roufines. 

We  follow  fhe  philosophy  of  our  previous  implemenfafions  [3],  and  choose  a  simpler  way,  based  on 
freafing  fhe  infinife  elemenfs  as  an  (implicif)  form  of  an  Absorbing  Boundary  Condifion  (ABC).  Insfead  of 
seffing  up  a  permanenf  dafa  sfrucfure,  we  generafe  a  temporary  data  every  fime,  fhe  problem  on  fhe  EE/IE 
mesh  is  solved.  The  resulfing  infinife  elemenf  d.o.f.  are  fhen  stored  in  fhe  femporary  dafa  sfrucfure  arrays, 
and  can  be  used  for  visualizafion  or  posf -processing.  The  momenf,  however,  any  finife  elemenf  adjacenf  fo 
fhe  fruncafing  surface  is  modified,  fhe  confenf  of  fhe  IE  dafa  sfrucfure  arrays  becomes  meaningless. 


The  femporary  IE  dafa  sfrucfure  arrays,  stored  in  module/infinite- elements,  include  fhree  arrays.  Arrays 


lELEMS  and  lENODS  store  new  objeets: 


•  type  infinite  element, 

•  type  infinite  element  node. 

The  attributes  of  an  infinite  element  inelude, 

•  the  element  nodes  listed  in  the  order:  middle  node,  left  radial  mid-edge  node,  right  radial  mid-edge 
node,  the  mid-edge  node  on  the  truneating  eirele,  left  vertex,  and  right  vertex  nodes, 

•  the  orientation  of  the  mid-edge  node  on  the  truneating  eirele. 

The  attributes  of  an  infinite  element  mode  are  similar  to  those  of  the  standard  nodes.  The  inelude  node  type 
(medg  or  mdlq),  order  of  the  node,  and  the  eorresponding  degrees  of  freedom  zdofs. 

The  third,  integer  array,  ELTOIE,  is  dimensioned  with  MAXNODS  -  the  maximum  number  of  FE  higher 
order  nodes.  The  array  is  populated  with  zeros  exeept  for  entries  eorresponding  to  middle  nodes  of  elements 
adjaeent  to  the  truneating  surfaee  whieh  store  the  adjaeent  infinite  element  numbers. 

The  edges  loeated  on  the  truneating  eirele  are  flagged  in  fhe  inpuf  file  wifh  (boundary  eondifion)  BC 
fiag  =  9.  The  infinite  elemenfs  are  generated  fhen  in  fhe  order  in  whieh  fhe  edges  wifh  fhe  BC  flag  =  9,  are 
eneounfered,  following  fhe  natural  order  of  elements.  In  ofher  words,  in  presenee  of  infinife  elemenfs,  fhe 
nafural  order  of  elemenfs  is  modified,  by  “inserfing”  infinife  elemenfs  after  fhe  adjaeenf  finife  elemenfs. 

Finally,  in  presenee  of  infinife  elemenfs,  fhe  inferfaee  wifh  fhe  fronfal  solver  musf  be  modified. 

Please  visif  fhe  souree  eode  for  defails. 

4.2  Calculation  of  infinite  element  stiffness  matrix.  The  infinite  element  routine 

As  diseussed  in  fhe  previous  seefion,  fhe  infinife  elemenfs  shape  funefions  are  tensor  produefs  of  FE  shape 
funefions  and  IE  radial  shape  funefions.  We  begin  wifh  a  paramefrizafion  for  fhe  edge  adjaeenf  fo  fhe 
fruneafing  eirele, 

a;e(6),6  e  (0,1).  (4.11) 

In  praefiee,  we  replaee  fhe  exaef  eirele  wifh  ifs  isoparamefrie  FE  approximafion, 

p-i 

®e(6)  =  *1(1  -6)  +  X]®3,fA3j(6)  •  (4-12) 

1=1 

Here  X3y(^i)  denote  fhe  mid-edge  node  (ID  master  elemenf)  shape  funefions  of  order  p,  a:i,a:2  are  fhe 
verfex  nodes  eoordinafes,  and  denofe  fhe  mid-edge  node  geomefry  d.o.f. 

The  mapping  from  fhe  infinite  master  element  oeeupying  fhe  sfandard  unif  square,  ean  fhen  be  deseribed 
as, 

*  =  *(6,6)  =  *e(6)^  • 

6 
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Notice  that  the  edge  ^2  =  0  corresponds  to  the  infinity.  The  IE  trial  shape  functions  are  defined  fhen  as, 

(4.13) 


Ak{x)  =  , 


where  fhe  radial  shape  funclions  are  defined  as  fhe  infegrals  of  fhe  Jacobi  polynomials  (fransformed 

fo  interval  (0,1)), 

'1  j  =  l 


= 


£ 

fi 


J  >  1  • 


(4.14) 


The  infinife  elemenf  fesf  shape  funclions  are  idenlical  wilh  fhe  Irial  shape  funclions  for  fhe  Bubnov-Galerkin 
version,  and  for  fhe  Pelrov-Galerkin  scheme  Ihey  are  given  by. 


(/>ika(x)  =  Xi(^l)  ^^4(^2), 


(4.15) 


where  fhe  radial  lesl  funclions  are  computed  according  lo, 

'  fo^  -  1)  dt 


'fpjiO  =  5 


,  /oPf’?(2f-l)df 


J  =  1 


J  >  1  • 


(4.16) 


The  compulation  of  Ihe  infinite  elemenl  stiffness  malrix  is  reduced  lo  evaluation  of  four  “one-dimensional” 


malrices 


{'^'Ipika'^fpji  -  k^-tpikafpji}  dx  =  a  ^  SjMli  +  a 


K 


(4.17) 


Here  S},,  5?  denote  “ID  curved  finite  elemenl”  malrices. 


cl 


c2 


x’MWoi^y 


(4.18) 


where. 


ds  dxg 

d^  d^ 

The  aclual  “infinite  elemenl”  malrices  are  given  by. 


Rh  =  -tp'k^i  dx  -  J 


Ki  = 


'fpkfl^i  d^  + 


—ika 


f 

Jo 


uo 

i’k(pi  d^ 


[  i'fPk^i 
Jo 


-'ipkfpi  )  dx  -  -ipkio)^^) 


(4.19) 


'Notice  that  integer  k  is  ‘bverloaded”. 


10 


Both  ID  stiffness  matrices  are  evaluated  using  standard  Gaussian  integration.  As  the  order  of  the  infinite 
elements  may  vary  from  an  element  to  element,  the  order  of  mid-edge  radial  nodes  is  determined  using  the 
minimum  rule.  The  element  stiffness  matrix  is  computed  first  for  the  element  of  a  “complete”  order,  with 
the  actual  element  matrices  extracted  from  the  complete  stiffness  matrix. 

Order  of  infinite  element 

The  infinite  elements  in  the  initial  mesh  are  assumed  to  be  isotropic,  i.e.  order  N  in  the  radial  direction  is  set 
to  the  corresponding  order  p  of  the  edge  on  the  truncating  circle.  The  BC  flag,  stored  for  fhe  corresponding 
mid-edge  nodes  is  sef  to  900  -|-  p. 

Typically,  we  begin  wifh  elemenfs  of  second  order.  During  fhe  /ip-refinemenls,  edges  on  fhe  fruncafing 
circle  gel  p-  or  /i-refined.  Every  lime,  fhe  edge  is  p-refined,  ils  IE  radial  order  is  updaled,  i.e.  fhe  BC  flag 
is  changed  from  900  -f  p  fo  900  -|-  p  -|-  1.  We  also  updafe  fhe  IE  order  when  fhe  edge  is  /r-refined  as  well, 
i.e.  fhe  sons  of  fhe  edge  are  assigned  fhe  modified  order  BC  flag  900  -|-  p  -|-  1.  Therefore,  in  presence 
of  /i-refinemenls,  we  encounter  infinite  elemenfs  wifh  radial  order  greater  lhan  fhe  EE  order.  This  rellecls 
fhe  philosophy  fhal  any  improvemenl  in  fhe  approximation  properties  on  fhe  fruncafing  circle,  should  be 
accompanied  wifh  fhe  corresponding  improvemenl  in  fhe  radial  direclion  as  well. 

In  fhe  presented  experimenls,  fhe  IE  radial  order  has  been  reslricled  lo  W  <  9. 


5  Calculation  of  RCS 

5.1  Direct  evaluation  using  the  IE  solution 

The  simplesl  way  lo  evaluate  RCS  is  simply  by  using  fhe  shape  functions  (4.13).  More  precisely,  if  fhe 
direclion  x  inlersecls  wifh  a  finile  elemenf  edge  (4.12)  on  fhe  fruncafing  circle,  and  is  fhe  value  of  fhe 
corresponding  parameler,  i.e., 

ke(6)r'"’ 

fhe  corresponding  formula  is  oblained  by  subsliluling  in  (4.13)  ^2  =  0, 

I  l®e(6)l  •  (5-20) 

k 

Here,  fhe  summation  extends  over  all  IE  shape  functions  (d.o.f.)  k,  {ik^jk)  denote  the  corresponding  indices 
for  the  ID  EE  shape  functions  Xi  and  radial  trial  functions  -ipj  constituting  the  IE  shape  function,  and  Uk 
denote  the  IE  d.o.f. 

If,  for  instance,  the  IE  formulation  is  stable  in  some  suitable  norm,  formula  (5.20)  should  give  convergent 
results. 
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5.2  Evaluation  through  postprocessing 


We  shall  limit  our  discussion  to  the  2D  case.  Any  function  u  that  satisfies  the  Helmholtz  equation  in  a  domain 
exterior  to  a  closed  boundary  F,  along  with  the  Sommerfeld  radiation  condition,  satisfies  automatically  the 
Helmholtz  representation  formula, 

f  f  Ou  ') 

u(x)  =  J^<l-  —  (y)^(x,y) +  u(y)  —  (x,y)j  ds{y)  .  (5.21) 

Here  n  =  n{y)  is  the  unit  normal  directed  to  infinity,  and  $  denotes  the  fundamental  solution  to  the 
Helmholtz  operator, 

^{x,y)  =  ^{k\x  -  y\)  , 


where. 


^{kr)  =  — ^ 


(2) 


{kr) 


4? 


(5.22) 


(2) 

Here  Hq  ’  denotes  the  Hankel  function  of  the  second  type,  of  order  0.  Notice  that  changing  the  ansatz  in 


time  to  e  requires  switching  to  the  Hankel  function  of  the  first  type. 


Evaluating  the  normal  derivative  of  the  fundamental  solution, 

0^  X  —  7J 

—  (a:,  y)  =  k  ^'{k\x  -  y\)  n{y)  ■  - — — r 
(jTZ  X  y 

' - V - - 

cos  fiy 


we  get. 

Using  the  far  field  approximation. 


y\)  +  k  u{y)^'{k\x  - 


(5.23) 


x-y\  =  \x\-y  -x  , 


and  the  asymptotic  formula  for  the  Hankel  function  and  its  derivative  [1], 


^{kr)  =  i 


^'{kr) 


^—ikr 


we  obtain. 


lim  |a:|2|u(a:)| 

\X  I  — ^oo 


1 

i' 


2 

ttA: 


-|^(y)  +  iku{y)n  ■  x 
IT  I  dn 


}e< 


k(y-x) 


ds{y) 


(5.24) 


The  integration  can  take  place  over  any  contour  surrounding  the  scatterer.  In  context  of  the  IE  discretization, 
it  is  natural  to  integrate  over  the  element  edges  adjacent  to  the  truncating  circle. 


6  Numerical  experiments 

In  all  reported  experiments,  the  error  is  computed  in  the  -semi-norm,  integrated  (only)  over  the  compu¬ 
tational  domain.  The  error  is  reported  in  percent  of  the  (semi)norm  of  the  solution,  defined  over  the  same 
domain.  The  exact  solution  for  the  cylinder  problem  is  computed  using  the  Mie  series.  Eor  the  wedge 
problem,  the  unknown  solution  is  replaced  with  the  solution  on  the  /ip-refined  mesh,  see  [6]. 
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Scattering  of  a  plane  wave  on  the  rigid  cylinder.  Verification  of  the  code 

We  begin  with  the  standard  eomparison  with  the  exaet  solution  for  the  problem  of  seattering  of  a  plane  wave 
on  a  unit  ey Under,  see  e.g.  [9].  We  set  the  wave  number  to  k  =  tt,  and  seleet  the  eomputational  domain 
in  the  form  of  a  ring  with  outer  (the  truneating  eirele)  radius  a  =  3.  Figure  2  displays  eonvergenee  history 
for  uniform  and  /ip-adaptive  rehnements.  As  expeeted,  the  uniform  p-rehnements  deliver  exponential  eon¬ 
vergenee,  with  the  adaptive  refinements  delivering  slightly  smaller  error  but  the  same  rates.  Figure  3  shows 
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Figure  2:  Seattering  of  a  plane  wave  on  a  eylinder.  Convergenee  history  for  uniform  p-  and  adaptive  hp 
refinements 

the  optimal  hp  mesh,  eorresponding  to  the  error  of  0.025  pereent.  Revealing  are  plots  shown  in  Figures 
4  and  5  presenting  eontour  lines  of  the  real  (range:  -0.62958E-02  0.57420E-02  )  and  imaginary  (range:  - 
0.47019E-02  0.45273E-02)  parts  of  the  error  funetion,  for  the  uniform  mesh  of  quartie  elements.  The  graph 
display  a  portion  of  the  infinite  elements  eorresponding  to  ^  <  ^2  <  1-  The  solution  in  the  IE  domain  is 
aetually  better  than  in  the  finite  element  domain  whieh  seems  to  be  eonsistent  with  the  faet  that  the  method 
seems  to  be  stable  in  norm  after  taking  out  the  Jacobian  (see  the  diseussion  in  Seetion  3).  Einally,  for 
the  same  mesh  of  quartie  elements  (1.5  pereent  error),  Eigure  6  presents  the  monostatie  RCS  eorresponding 
to  range  from  180  to  0  degrees  (left  to  right),  displayed  in  dB  (20  log  of  the  aetual  value).  Both  methods 
for  eomputing  RCS  yield  identieal  results  (first  three  signiheant  digits  aeeurate).  The  derivatives  in  formula 
(5.24)  have  been  averaged  over  the  adjaeent  Unite  and  inUnite  elements.  With  the  derivatives  evaluated  using 
only  the  eontributions  from  either  Unite  or  inUnite  elements,  the  eorresponding  value  of  RCS  differs  in  the 
third  digit.  In  this  sense  the  Urst  approaeh  seems  to  deliver  better  results  whieh  is  quite  remarkable.  The 
slight  variations  of  the  RCS  for  the  eylinder  reUeet  the  imperfeet  approximation  of  the  geometry. 
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Figure  3:  Scattering  of  a  plane  wave  on  a  cylinder.  Optimal  hp  mesh  corresponding  to  0.025  percent  error 

Comparison  of  Bubnov-  and  Petrov-Galerkin  formulations 

Fig.  7  displays  the  convergence  history  for  the  automatic  /ip-adaptivity,  and  the  two  versions  of  the  infinite 
element  test  functions.  The  Petrov-Galerkin  version  delivers  slightly  better  results  for  coarse  meshes  but  the 
differences  are  insignificant.  Both  methods  seem  to  display  the  same  character  of  convergence.  In  particular, 
the  RCS  computations  discussed  above  yield  identical  results  (first  three  digits  coincide).  Same  conclusions 
apply  to  the  wedge  problem  discussed  next. 

Scattering  of  a  plane  wave  on  a  wedge 

The  second  and  final  example  deals  wifhe  fhe  resolufion  of  singularifies,  and  ifs  impacf  on  RCS.  We  have 
divided  fhe  cylinder  from  fhe  previous  example  info  four  equal  quadranfs  and  kepf  jusf  one  of  fhem.  We  sef 
fhe  wave  number  k  =  tv,  i.e.  fhe  disfance  befween  fhe  fruncafing  boundary  and  fhe  objecf  is  equal  fo  one 
wavelengfh.  Jusf  as  an  example,  Figures  8  and  9  presenf  fhe  optimal  hp  mesh  wifh  zooms  on  one  of  fhe 
comers  for  fhe  rafher  academic  error  level  of  0.1  percenf.  Fig.  10  presenfs  fhe  convergence  hisfories  for 
fhe  problem  using  bofh  Bubnov-  and  Pefrov-Galerkin  formulations.  Similarly  fo  fhe  cylinder  problem,  fhe 
differences  are  insignificanl  allhough,  again,  fhe  Pelrov-  Galerkin  melhod  seems  fo  be  doing  a  lillle  heller. 
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Figure  4:  Scattering  of  a  plane  wave  on  a  cylinder.  Real  part  of  the  error  function  for  a  uniform  mesh  of 
quartic  elements 

Does  the  resolution  of  singularities  matter  ? 

We  come  to  the  final  experiment  reflecting  the  impact  of  adaptivity  on  evaluation  of  RCS.  Figure  1 1  presents 
RCS  for  the  wedge  problem  evaluated  using  a  uniform  mesh  of  quartic  elements  and  /ip-adaptive  meshes. 
The  choice  of  the  uniform  meshes  reflects  the  usual  practice  in  selecting  a  mesh  that  reproduces  the  wave 
form  of  the  solution  (two  quartic  elements  per  wavelength)  and  delivers  an  error  in  the  range  of  3-4  per¬ 
cent.  The  hp  meshes  were  obtained  by  requesting  a  one  percent  error  level,  at  which  several  levels  of 
/ip-refinements  resolve  the  structure  of  the  singularities  in  the  solution.  For  each  direction  of  the  incoming 
wave  (0  =  180, 179, . . .  ,  0,  left  to  right),  the  /ip-adaptive  algorithm  was  mn  starting  with  the  optimal  mesh 
for  the  previous  direction,  with  the  optimization  procedure  restarted  from  the  initial  mesh  every  10  degrees. 
The  result  is  discouraging  for  those  “selling”  adaptivity.  Except  for  a  slight  shift  in  the  RCS  level,  the  results 
are  practically  identical.  Resolution  of  the  singularity  seems  to  be  insignificant. 

On  the  positive  (for  adaptivity)  side,  the  same  error  level  of  roughly  3.5  percent,  and  the  quality  of 
the  corresponding  RCS,  can  be  obtained  with  meshes  roughly  half  the  size  of  the  uniform  meshes.  Figure 
12  presents  an  example  of  such  a  mesh  for  the  direction  of  90  degrees.  Again,  the  fact  that  the  algorithm 
selected  only  p-refinements,  indicates  that,  at  this  error  level,  the  resolution  of  singularities  does  not  matter. 
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Figure  5:  Scattering  of  a  plane  wave  on  a  cylinder.  Imaginary  part  of  the  error  function  for  a  uniform  mesh 
of  quartic  elements 

Evaluation  of  RCS 

Both  methods  of  evaluating  RCS  seem  to  deliver  practically  the  same  results  and  converge  with  the  same 
rates.  The  direct  method  seems  to  be  doing  a  bit  better  as  illustrated  in  Fig.  13. 

7  Conclusions 

We  summarize  the  main  conclusions  of  our  brief  considerations  and  numerical  experiments. 

1 .  Interpreting  the  calculation  of  the  sesquilinear  form  in  the  CPV  sense,  allows  for  a  symmetric  choice 
of  trial  and  test  shape  functions.  The  choice  does  not  seem  to  have  any  impact  on  the  convergence 
properties.  Will  the  symmetric  choice  be  easier  to  analyze  in  terms  of  convergence  proofs? 

2.  Both  formulation  allow  for  a  direct  evaluation  of  RCS  from  the  IE  representation  of  the  solution.  The 
method  seems  to  converge  uniformly  in  the  whole  exterior  domain. 

3.  The  use  of  locally  variable  order  of  IE  approximation  is  possible.  The  infinite  element  can  be  used  in 
conjunction  with  the  automatic  h.p-adaptivity. 

4.  Resolution  of  singularities  seem  to  have  little  impact  on  the  RCS  calculations,  although  adaptivity 
(ideally  goal-driven)  allows  for  essential  savings  in  terms  of  problem  size. 
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2.74  RCS 
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Figure  6:  Scattering  of  a  plane  wave  on  a  cylinder.  RCS  for  the  uniform  mesh  of  quartic  elements  (1.5 
percent  error) 

5.  The  most  important  difference  between  the  adaptive  and  non-adaptive  techniques  comes  with  the  fact 
that  the  adapted  mesh  comes  with  an  error  estimate  guaranteeing  the  quality  of  the  delivered  numerical 
result. 

The  hp  code  used  in  the  presented  experiments  can  be  downloaded  from 

http://www.ticam.utexas.edu/%7Eleszek/hp-intro.html 
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Figure  8:  Scattering  of  a  plane  wave  on  a  wedge.  Optimal  hp  mesh  for  0.1  percent  error,  with  a  10  times 
zoom  on  the  lower  comer 


Figure  9:  Scattering  of  a  plane  wave  on  a  wedge.  Optimal  hp  mesh  for  0.1  percent  error.  Zooms  on  the 
lower  corner  with  100  and  1000  magnifications 
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Figure  10:  Scattering  of  a  plane  wave  on  a  wedge.  Convergence  history  for  adaptive  hp  refinements  using 
Bubnov-  and  Petrov-Galerkin  IE  formulations 


Figure  11:  Scattering  of  a  plane  wave  on  a  wedge.  RCS  for  the  uniform  mesh  of  quartic  elements  (3-4 
percent  error  range  level)  and  /ip-adaptive  mesh  (1  percent  error) 
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Figure  12:  Scattering  of  a  plane  wave  on  a  wedge.  Optimal  hp  mesh  for  the  direction  of  90  degrees 


Figure  13:  Scattering  of  a  plane  wave  on  a  cylinder.  Convergence  of  the  RCS  evaluated  using  the  two 
methods 
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